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Abstract 

For numerical simulations of highly relativistic and transversely accelerated charged particles including 
radiation fast algorithms are needed. While the radiation in particle accelerators has wavelengths in the 
order of 100 urn the computational domain has dimensions roughly 5 orders of magnitude larger resulting in 
very large mesh sizes. The particles are confined to a small area of this domain only. To resolve the smallest 
scales close to the particles subgrids are envisioned. For reasons of stability the alternating direction implicit 
(ADI) scheme by D. N. Smithe et al. (J. Comput. Phys. 228 (2009) pp.7289-7299) for Maxwell equations 
has been adopted. At the boundary of the domain absorbing boundary conditions have to be employed to 
prevent reflection of the radiation. In this paper we show how the divergence preserving ADI scheme has to 
be formulated in perfectly matched layers (PML) and compare the performance in several scenarios. 
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1. Introduction 

For simulations of electromagnetic fields an implementation using finite differences suggests itself through 
its simplicity and low computational expense. The explicit Yee scheme pQ is very well studied and widely 
adopted. Although its explicit formulation is to its advantage when it comes to computing time per time step 
it is also its weakness in terms of numerical stability. The Courant-Friedrich-Lcwy (CFL) condition limits 
the ratio of the length of a time step and the width of a mesh cell. If high spatial resolution is needed only 
in a small fraction of the computational domain, then this constraint causes computational overhead, since 
the smallest cell determines the largest time step to be used. To overcome this constraint on the time step 
an alternating direction implicit (ADI) scheme can be adopted [5J [3] . In this scheme the time step is split in 
two substeps. For every component of the electromagnetic field the partial derivative, originating from the 
curl operator, with respect to one direction is updated alternately implicitly and explicitly while this choice 
is reversed for the other partial derivative. In 2D half of the partial derivatives vanish, thus simplifying 
the systems of equations. The implicit part of the equations for the electric and magnetic components are 
coupled in such a way that multiple independent systems of equations in only one dimension emerge. For the 
simplest case of a homogeneous linear medium and a perfect electric conductor surrounding it, the resulting 
matrices are symmetric tridiagonal with dominant diagonals. Hence the resulting system is stable and can 
be solved fast. This amounts to a Cholesky decomposition at the beginning of a simulation and forward and 
backward insertion in every time step and spatial direction. 
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In simulations of charged particles that include the induced electromagnetic fields it is important to 
calculate the current and the fields self-consistently. To achieve self-consistency the interpolation of the 
electromagnetic field a the position of a particle and the interpolation of the charge current at the vertices 
of the mesh have to preserve the charge conservation, V • J + J^p = 0, in a discrete sense. Furthermore, 
the generated fields have to satisfy Maxwell's equations. The latter gives rise to the necessity of a preserved 
divergence of the fields (Gauss's law). Smithe et al. [I] showed that the original ADI scheme for Maxwell's 
equations introduced in [5] does not satisfy this requirement. After a thorough study of all second-order 
ADI schemes they derived the only ADI scheme for Maxwell's equations that does preserve the divergence. 

For simulations of open domains absorbing boundary conditions with low reflection are used to confine 
the computational domain. Perfectly matched layers [5] accomplish absorption by introducing an electric 
and a magnetic resistivity in layers at the boundary of the domain. The Maxwell's equations in these layers 
then become 

Q 

e -E + aV = V x H, 
at 

d 

Mo^H + cr*H = -V x E, 
at 

where a and a* are the electric and magnetic conductivities, respectively. The PML can be adapted to the 
ADI scheme [BJ. To that end the conductivity term is split and added equally to both the new and the old 
time step. For large time steps this choice is not optimal as was shown in [7]. The reflections can be reduced 
if the conductivity term is not split in equal parts. Instead it is much more advantageous to add it in one 
substep and in one direction fully to the new time step while in the other direction it is added to the old 
time step. In the next substep this choice is then reversed. 

In the following section we will shortly introduce the original and the divergence preserving ADI scheme 
in the formulation of [I]. We restrict our presentation to the 2D case. Then we will show how perfectly 
matched layers are formulated using the divergence preserving scheme. To confirm the validity of the 
resulting formulation we conduct similar numerical experiments as in [7] and compare the results with the 
ones produced with the Yee scheme and the original ADI scheme in Section [3] 



2. Alternating Direction Implicit Schemes 

Using the four operators that were elaborated in [1] the original ADI scheme for Maxwell's equations 
reads 

1 + C -^v) 1 ■ (l C -^m) ■ V-+ 1 = (l - C ^V) 1 ■ (l + C ^M ) ■ V™ + c At 



where in 2D 







E y 
ZoH" J 


I- 




f Zq J x 


-1/2 


Z Q Jy 


-1/2 







Z V y H z /Ay 



v ■ v n = i v 



V f y E2/Ay 

S"" 1 2 = - I Zn ./r 1 2 I • M>V n = -\ Z Q V b x H z /Ax 



HKIA 



x 



Here, Zq is the impedance and c the speed of light in vacuum, J the current density and T>[, T) h x the forward 
and backward difference operators in x-direction, respectively. Their action on E y is described by 

, = E y . l+ ij - Ey. i} j and 
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The operators T>' y and D b are the corresponding operators in y-direction. In this notation the curl operator 
is(M + P). 

The divergence preserving ADI scheme is built up using the same set of operators but applied in a 
different order: 
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Remark: In 3D the vectors V, S and the operators A4, V have to be adapted: 
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3. Perfectly Matched Layers 

The absorption of electromagnetic waves in PML is accomplished by introducing layers of artificial media. 
Each medium is described by an electric and a magnetic conductivity. The electric conductivity, a, causes 
an exponential reduction of the amplitudes of incoming waves as e~ a m z ° x . The magnetic conductivity, 
a*, on the other hand is chosen such that the impedance of the medium matches the one of vacuum. To 
reduce the numerical reflections caused by abrupt changes in the conductivities, media with increasing 
conductivities towards the boundary of the domain are applied. The boundary of the domain is then 
terminated by a perfect electric conductor (PEC) . The total theoretical reflection of a plane wave is described 

by R((f) = e -2^oco S¥ >/ *(r)dr 

if is the angle between the normal of the boundary and the direction 
of propagation of the incident wave and 5 is the width of the PML region [5]. 

In 2D the electric conductivity is split in two components, a x and a y , which are associated with the y- 
and x-component of the electric field respectively. Also the magnetic field is split in two subcomponents, H zx 
and H zy (in 3D we have an additional four components of the magnetic field H xy , H xz , H yx , H yz and two 
components of the electric field E zx , E zy while the two existing components are split in four subcomponents 
E xy , E xz , E yx , E yz ). Finally, the magnetic conductivity is split in components, a* and <r*, associated with 
H zx and H zy , respectively. 

In our approach we add the conductivities a x , a y , a* and a* to the four operators in ([T]). The divergence 
preserving ADI scheme including PML then reads 

p-yl^T^'^^^yl^-f^ 1 ■v" + C Ats»+v 2 , (2) 
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PML 

-> E x t E y © H z • H zx & H : 



Figure 1: Example of a domain terminated by five layers of perfectly matched media. The increase of cr y and cr* towards the 
boundary of the domain is visualized by the shading. 
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and, for r = {1, 2, 3, 4}, 
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With A|* and A|' J one can control whether the conductivity terms are added fully to quantities at the new 
time step (+1) or to the quantities at the old time step (0) or a linear combination thereof. 
The first substep of our new scheme inside the PML reads 
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The second substep is then given by 
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As one can see from the inverted operators in ^ and Q the systems that have to be solved are only in one 
dimension and the resulting matrices are tridiagonal. This is due to the fact that V couples E x and H zy 
only and the operator M. couples E y and H zx . In 3D a one-dimensional system has to be solved for every 
component of the electric field. 

In [7] it is shown that for the non-divergence-preserving scheme one gets the best result if in one substep 
the conductivity is added fully to the quantities at the new time step of the pair (E x , H zy ) while it is added 
fully to the quantities at the old time step of the other pair, (E y , H zx ). In the second substep it is the other 
way round, for (E x , H zy ) it is added fully explicitly and to (E y , H zx ) it is added fully implicitly. For the 
divergence preserving scheme we find the combination 

A x 1 = \' 1 = 0; 

\e, 2 _ \h, 2 _ 1 
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to yield the lowest reflections. This choice for A?' T adds the conductivity term to the new time step if the 
operator in which it appears is applied implicitly and adds it to the old time step if the operator is applied 
explicitly. In 3D the set of a's is augmented by the a's of the additional components of the electromagnetic 
field. The rule how to choose the A's is however not altered. 
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4. Numerical Experiments 

To validate the performance of PML in the divergence preserving ADI scheme a test domain with 21 
cells in each direction is excited by a source at the center, see Fig. [2] The domain is enclosed by 10 perfectly 
matched layers on each boundary with 



where r is the distance to the interface between vacuum and PML. a rn is chosen such that the theoretical 
reflection at an incident angle of is e -16 [H Ch. 7.6]. The source is a sine wave with frequency 3.175 GHz 
modulated by a differentiated Gaussian pulse with half bandwidth of 3.172 GHz [7], 

H z (x , j/o ) = H z (x , j/o ) - — s-^ sin(w • t) e~ ' ^ , 

K 2 -l 

with u = 2ir ■ 3.175 • 10 9 , t = 0, t = ^ and n e — = 0.01. The simulation is run for 1.5 ns or - 4.75 
periods of the central frequency. 
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Figure 2: Setup of the numerical experiment with a soft source in the center of the domain, a probe near the interface to the 
PML and 10 Layers of PML on the boundary of the domain. 



A probe to sample the magnetic field is placed 10 cells away from the source, next to the interface between 
vacuum and PML. As a reference signal the magnetic field at a single point in a larger domain is recorded. 
This reference domain has 411 cells in each direction and is terminated by PEC boundary conditions but 
has no PML. The source is also placed at the center of the domain and the location of the probe is 10 cells 
away. The size of the reference domain ensures that the backscattered waves do not reach the probe during 
the simulation. The relative difference of the readings of the two setups, 

\H z (t)-H z . Icl (t)\ \ 
max\H z . rcf (i)| J ' 

versus time is plotted. The length of the time step in the simulation relative to the maximal time step for 
the Yee scheme is characterized by the Courant number \ — Ai/Ai max where At max = Ax / (y/2v p ) and v p 
is the phase velocity. 

4-1. Small Time Steps 

For a Courant number x = 0.5 the reflections for the divergence preserving ADI scheme and the Yee 
scheme are depicted in Fig. [3] on the left. On the right side data of a simulation with \ — 1, i.e. the largest 
time step feasible with the Yee scheme, are depicted. It is clearly seen that for small time steps the PML for 
the divergence preserving ADI scheme performs equally well as the Yee scheme. For x = 1 the Yee scheme 
performs slightly better but the results of both schemes are still below —90 dB. 

4-2. Large Time Steps 

For x = 6 the divergence preserving scheme is compared with the original ADI scheme. The Yee scheme 
cannot be used at such large time steps since it becomes unstable for x > 1. As one can observe in Fig. [4] 
(left) the two schemes produce the same result up to relative differences in the order of 10 -10 . The absolute 
difference, Aadi. dp-Adi = Radi — -Rdp-Adi, is added to the plot. If one compares the reflections of the ADI 
scheme with x — 6 in Fig. [2] (left) with the Yee scheme with x — 1 m Fig. [3] a considerable degradation of 
the performance of the PML can be observed. Nevertheless, the reflection is still below —60 dB, which is 
sufficient in many situations. In Fig. 4 (right) the divergence preserving ADI scheme with A^' 1 = A^' 3 = 0.5 
(conductivities added half implicitly, naif explicitly) is compared with the same scheme but the proposed 
A's for x = 0.5 and x = 6. As was already observed in [7] an improper choice of the weights A^ r in ([2| 
leads to a comparable degradation as increasing the step size by a factor of 12. 
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Figure 3: Comparing the reflection at layers of perfectly matched media of the Divergence Preserving ADI and Yee scheme for 
X = 0.5 (left) and \ = 1 (right). To reduce the noise the data was smoothed using a rolling average. For the case \ = 0.5 a 
length of 10 data points is used and only every 10 th point is drawn. For the case \ = 1 om y half of this length is used for the 
window and every 5 th point is drawn 
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Figure 4: Comparing the reflection at layers of perfectly matched media of the Divergence Preserving ADI and the original 
ADI scheme for x = 6 on the left. On the right side the reflections of four runs with the divergence preserving ADI scheme are 
compared. Two of these runs were carried out with the proposed A's (one-sided) while for the other two runs the conductivities 
were added half implicitly and half explicitly (equally). For both choices of the A's the simulation was carried out with two 
different lengths of the time step, i.e., \ = 0.5 and x = 6. 
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5. Conclusions 

It has been shown how perfectly matched layers have to be formulated within the divergence preserving 
ADI scheme. Furthermore their performance has been validated by comparing the reflections for small time 
steps with results obtained with the Yee scheme. This comparison gives very good results. 

For large time steps beyond the CFL limit the implementation has been compared with results using a 
non-divergence preserving ADI scheme. The divergence preserving scheme gives equally low reflections as 
the non-preserving scheme. 

Perfectly matched layers combined with a divergence preserving ADI scheme allow to perform simulations 
of charged particles in (partially) open domains. With this implicit scheme, the time step in these simulations 
is not bounded by the CFL condition. In our implementation the computing time for one step on a single 
core was roughly three times as large for the ADI scheme as for the Yee scheme. So, if the time step in the 
ADI scheme is chosen three times larger than in the Yee scheme then the ADI scheme gives faster results 
with sufficient accuracy. 
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